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We have studied the quantum phase transition between the antiferromagnetic and spin liquid 
phase for the two dimensional anisotropic Kondo necklace model. The bond operator formalism has 
been implemented to transform the spin Hamiltonian to a bosonic one. We have used the Green's 
function approach including a hard core repulsion to find the low energy excitation spectrum of 
the model. The bosonic excitations become gapless at the quantum critical point where the phase 
, transition from the Kondo singlet state to long range antiferromagnetic order takes place. We have 

^ . studied the effect of both inter-site (5) and local (A) anisotropies on the critical point and on the 

' critical exponent of the excitation gap in the paramagnetic phase. We have also compared our 

results with previous bond operator mean field calculations. 
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The description of quantum phase transition between phases with spontaneously broken symmetry and disordered 

S phases is a novel topic in condensed matter physics^. Macroscopic strongly correlated electron systems at low 
temperature (and as a function of magnetic field, hydrostatic or chemical pressure) show a wide range of interest- 
, ing phenomena, such as quantum criticality and associated non-Fermi liquid (NFL) behaviour, magnetism, Kondo 
insulating behaviour and superconductivity^. In the single-impurity case the Kondo problem describes the antifer- 
romagnetic interaction (J) between the impurity spin and the free conduction electron spins. This gives rise to a new 
non-perturbative low-energy scale, the Kondo temperature Tk which dominates the low temperature anomalies in the 
thermodynamic and transport quantities 3 . Tk = De~ 1 ^ 2Jp ^ (D, p are conduction band width and density of states, 
H ■ respectively) has the meaning of a crossover temperature from uncoupled local spins for T 3> Tk to the strongly 
>■ coupled local spins, forming a singlet ground state with conduction electrons, for T <c Tk- In the Kondo lattice (KL) 
^sO , models an additional (perturbative) energy scale Trkky = J 2 p for the effective Ruderman-Kittel-Kasuya-Yosida 
• (RKKY) inter-site interactions of local spins appears. This model exhibits a quantum phase transition between the 
Kondo singlet phase and the magnetically ordered phase as function of the control parameter x = J p{Ep) as argued 
by DoniacbA The transition takes place at a quantum critical point (QCP) characterised by x c = J c p{Ep) where 
x c is of the order one. For x <C x c the effective interactions dominate and magnetic order appears. For x ^> x c the 
singlet formation dominates and a heavy Fermi liquid state is realized. This qualitative picture has been supported 
by numerical calculations within dynamical mean field theory (DMFT) and exact diagonalization methods^ for the 
Anderson lattice and Kondo lattice Hamiltonian respectively. However the vicinity of the quantum critical point and 
associated NFL behavior— requires a treatment within phenomenological effective models as developed in Refs. [TollTTI . 

The Kondo lattice model emerges from periodic Anderson model via a Schrieffer- Wolff transformation that elimi- 
nates the fluctuations of f-charge or f-orbital occupation^ . It is given by 
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H K L = t ( C L C 3,r + h.C.) + J±J2 TA ■ W 



The first part describes conduction electrons c| T with n.n. hopping t. The second part is the Kondo term where 
Tj and Si are conduction electron and localized spin respectively. This model still contains the charge fluctuations 
of conduction electrons expressed by the hopping term. It was shown by DoniacbJi that in ID it may be replaced 
by an xy-type inter-site exchange term. Thus the Kondo lattice model is replaced by a pure spin Hamiltonian, the 
Kondo- necklace model (KNM). In higher dimension this procedure cannot be justified strictly. However suppose we 
add a Coulomb repulsion U c between conduction electrons to the KL model (Eq. [1]) in the half filled case. Then in 
the limit U c /t — > oo charge fluctuations of conduction electrons are frozen out and the low energy physics is again 
described by a pure spin Hamiltonian. Strictly speaking this is only adequate for the Kondo insulator with a charge 
gap but one may expect that it is also useful to describe the low energy spin dynamics of metallic Kondo systems. 
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The generalized Kondo necklace model obtained in this wayi£ is given by 

H = J $>fr; + tV tJ + drfr?) + J ± £(rf Sf + rf Sf + ArfSf) , (2) 

(hi) <*) 

where the intersite exchange J is of order P /U c . In 2D which will be considered in the present work this is equivalent 
to a special case of an (anisotropic) bilayer-Heisenberg models where the inter-site bonds (J) are cut in one layer and 
J_l is the inter-layer coupling. Here both spins are 1/2 and the exchange coupling parameters are antiferromagnetic 
(J, Jj_ > 0). In the above Hamiltonian, r" represent the a-component of spin of the 'itinerant' electrons at site i and 
Sf is the a-component of localized spins at position i. 

We want to study the possible quantum phase transition of this model under rather general assumption of both 
anisotropies in the inter-site interaction (~ J) and on-site Kondo terms (~ Jj_) of Eq. They are characterized by 
a pair of parameters (6, A). The A anisotropy is always present in real Kondo compounds like Ce-based intermetallics 
due to the crystalline electric field (CEF) and 5 is caused by spin-orbit coupling of conduction electrons. We study 
the quantum phase transition from the paramagnetic (Kondo-singlet) side as function of the control parameter Jj_/ J, 
which gives the ratio of the inter-site to the on-site interaction strength, and as function of the anisotropy parameters 
(6, A). We have implemented the Green's function approach introduced to study the bilayer isotropic Heisenberg 
model (S, A) = (1, 1)^2. The effect of anisotropies on the quantum phase transition of KNM has recently been studied 
by a mean field approach both in the absence^ and presence^ of a magnetic field. However, using the more advanced 
Green's function method we will obtain more accurate values for the critical gap exponents v and for quantum critical 
point values (J±/J) c which differ both from the mean field values. 



II. BOSON OPERATOR REPRESENTATION OF THE MODEL HAMILTONIAN 

The bond operator representation introduced by Chubukov^ and Sachdev, et.ali 7 - is a useful approach to describe 
disordered phases. This representation can be considered as an analog of the usual Holstein-Primakov transformation 
for phases with broken spin rotational symmetry. In terms of singlet-triplet operators, the spin operators of the 
localized and conduction electrons are given by 
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(3) 



where (a,/3, 7) represent the (x,y,z) components and e is the totally antisymmetric tensor. The bond operators 



satisfy bosonic commutation relations 



u. 
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1. We will calculate the one 



particle boson Green's function using Feynman diagrams for triplet operators and find the excitation spectrum. Our 
calculations are for zero temperature. In order to ensure that the physical states are either singlets or triplets one has 
to impose the constraint s + ^ Q P a t a = 1 on every bond where s (singlet) and t a (triplet) are bond operators. 

In Refs. IT^fTll this has been implemented on a mean field level by introducing a chemical potential as Lagrange 
parameter. In this approach average amplitudes s = (si) and i = (U a ) are introduced and their self consistent solutions 
are found by minimizing the total ground state energy. Here t 7^ denotes a triplet condensed state with magnetic 
order. The chemical potential adjusts itself such that the averge constraint s 2 — 1 — P is approximately satisfied. In 
fact in the paramagnetic region (t = 0) it was foundM that s is only a few per cent below singlet saturation s = 1 even 
close to the QCP where triplet excitations become soft. The zero point energy of the latter contribute to the ground 
state energy. Since in the mean field approach the number of triplet bosons on a given bond is not constrained there 
are contributions from unphysical states in the ground state energy. 

In the present work we are therefore using a more advanced implementation of the local constraint which can 
be written as s]si = (1 — ^2 a t\ a ti a ). It may be satisfied if either sjsi = 1 and ^2 a t\ a ti a = or sjsj = and 
^2 a t\ a ti a = 1. To project out unphysical states on every bond with more than one excited triplet one has to require 
Aa^la = 0- This may be achieved by introducing an on-site repulsion U of triplet bosons^ which is then taken in 



the hard core limit U — > 00, see Eq. (|12p below. As starting point for nonintcracting triplets we use the unconstrained 
case with s — > s — 1 in Eq. ([3]). In the paramagnetic case which we consider here this is well justified by the mean field 
result mentioned above. In principle one might think of a combined approach keeping s as a variational parameter 
within the hard core boson approximation. We will discuss this further in Sect. rVTHl This hard core boson approach 
can be applied to any model, for which the excitations in the disordered phase are triplets above a strong coupling 
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singlet ground state. The Hamiltonian in Eq. ([2]) has three control parameters, J, 5 and A. Using the bond operator 
transformations in the Kondo-necklace model of Eq. @ , we obtain the effective Hamiltonian 

H = H 2 +H 3 + H i , (4) 

where H 2 is the one particle part of the Hamiltonian. It is composed of two terms 

H 2 = H J± +H 1 , (5) 

where the exchange term Hj ± is diagonal in terms of the bond operators and H± has pairing terms between boson 
triplets which results in non-conservation of the triplet bosons and the possible formation of a Bose-Einstein condensate 
of triplet bosons describing the magnetically ordered state. This term leads to a nonzero anomalous expectation value 
(tat a) ° r corresponding anomalous Green's function. In terms of bond operators, Hj ± and Hi are given by 



H 



J± [t\, x ti, x + tlj i<y ] + t{ z U iZ ) , (6) 



(id) ot=x,y (ij) 



The other parts of the Hamiltonian which describe triplet boson interactions are represented by 

J 

(id) 

+ S(U,, + tl z )(t), x t 3 ,y - *},„«,>)] + /i.e.) , (7) 



<i,j) 



The Hamiltonian can be written in terms of triplet Fourier components, t i a — Y] h a t kM e l { , leading to the 
quadratic form 

H 2 = A k J k J Ka + ^f(i,J-k, a + h.c). (9) 

k,a=x,y,z k,a—x,y^z 

The coefficients in the above equation are 

A k , z = J±+SJCk , A k:(Xiy) =^-(l + A) + JCk, 

Bk,z=SJ^k , B kt ( XtV ) = J£k, 

£/- = {cosk x + cosk y ]/2. (10) 

Also for H 3 we obtain 

H3 = iJ E £ fel (^x,k 1 ^y,k 2 ^ z M ~ *l,fci*l,fe2*3/. fe 3 + ty tkl t z j. 2 tx,k 3 

fcl ,fc 2 ,fc 3 = fei +/c 2 



Because f/3 and i?4 are of higher order in triplet operators they will lead to only small corrections in the spectrum. 
Therefore the effect of H 3 and H4 may be taken into account on a mean field level. 
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The dominant contribution to the renormalization of the spectrum comes from the constraint where only one of 
the triplet states can be excited on every site (the hard-core condition) f^tL = 0, which can be taken into account 
by introducing an infinite on-site repulsion between the bosons 

H V = UJ2 4i*l**/M*«. U —> °°- ( 12 ) 

Writing Hjj in terms of Fourier transforms of boson operators we obtain 

H U = U ]T t\ k+q tl k ,_ q t 0k d ak . (13) 

k,k' ,q,ct,(3 



III. GREEN'S FUNCTION FORMALISM IN THE BOSONIC TRIPLET GAS 

The second part of H 2 leads to the non-interacting normal Green's function and in addition to anomalous Green's 
function. Therefore, we introduce the single particle Green's function for the non-interacting Hamiltonian which help 
us to obtain the interacting {Hjj) Green's function for the triplet operators by using Dyson's equation. Implementing 
the Bogoliubov transformation tk, a = u kyCc t k ^ a + Vk >a F_ k we obtain w ka = A\ — B k for the excitation spectrum 
at the quadratic level (H2 only). The Bogoliubov coefficients are u k a (v k a ) = (— )| + ■ The non-interacting 
normal triplet Green's function is G^{k,t) — —i(T(t k:Cl (t)t k (0))) and the anomalous Green's function is G^(k,t) = 
-i(T(tl,a(¥-k, a (°)))- Together we have 

G n a (k,u) = ^— k > a - (14) 

W - U>k,a +IV U + U) kta - IT) 

G (k,uj) — ■ — ; — . (15) 

OJ - UJk.a + IT} U + ^fe,a - If] 

The interacting Green's functions are obtained from Dyson's equation for each Green's function (anomalous or normal). 
The perturbation expansion for the interacting Green's functions (for each polarization component of the triplet 
bosons) is written by 

G(k,u) = GP(k,w)(l -G°(fc,w)E(fc,w)) _1 . (16) 
The interacting Green's function (G(k,ui)) and the self-energy (S(fc, ui)) are 2x2 matrices 

B <*-> = ( oTA ) • = ( ) ■ M 

The matrix form of Green's function can be simply expressed by G(k,t) — (T($(fe, i)$t(fe, 0))), where t) = 

(tf(k,t) t(—k,t)) is the row vector. Inserting the elements of Eg . (|T4"j) , Eq. lfTB")) into Eq. lTBl) . the normal and 
anomalous interacting Green's function will be obtained by 



G n , a (k,uj) = 



[w + A k , a + S„,a(fc, -u)][uj - A k _ a - T, n>a (k,uj)] + (B k + E at0l (k,u>)) 



2 ' 



a ' a ' [cj + A k ^ a + £„, a (fc, -u>)][uj - A k ^ a - S„, a (fc,w)] + (B k + T, a . a (k,oj)) 2 ' 

The normal and anomalous self-energy are due to boson interactions #3 and Hjj and will be discussed in the next 
sections. Here we are interested to find the one particle excitations which are the poles of the normal triplet Green's 
function. The Green's function should be separated into the bosonic excitation contribution and incoherent back- 
ground (including collective modes) . Indeed, poles of the one particle Green's function of the triplet bosons result in 
low energy excitations of the Hamiltonian which vanish close to the critical point. To get the single particle excitation 
the self-energy is expanded for low energies leading to 

r , _ w + A k ^ a + E n , a (fc, 0) - ud u T; n . a (k, 0) 

"^'^ D + {B k . a + ^ a {k,W 

D = D!-D 2 , 

Di = [w + il|fc,a + Iln,a(A,0)-wS w E„ ) a(A,0)], 

D 2 = [lu- A Ka - 0) - wd u En, a {k,0)]. (19) 
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Splitting Eq. ([T9|) into partial fractions leads to the single particle (sp) parts 

G^ Q (fc,w) = - '— ! — — , (20) 

where the renormalizcd triplet spectrum and the renormalized single particle weight constants are given by 



7-1 _ 1 ( ^n,a v 

TT 2 / t/ 2 \_[ \1 | Zk y a[Ak y a + E ra! q(fc, 0)] . . 

u k,a\ V k,a) ~ rJn ' ^ ■ 

The renormalized weight constant is indeed the residue of the single particle pole in the Green's function. For the 
non-interacting system it is equal to one. 

IV. CALCULATION OF BOSON SELF-ENERGY DUE TO H v AND H 3 

Since the Hamiltonian Hu in Eq. (fT3]) is short ranged and U is large, the ladder diagram approach^ may be 
applied. This approach is suitable to solve Dyson's equation in order to get the boson Green's function. Formally 
this is quite similar to Ref. [l3j, however, technically more demanding due to the effect of anisotropies (<5, A). 

Now, we should impose the hard core repulsion due to the Hamiltonian Hu and obtain the interacting normal 
Green's function by Dyson's equation. Firstly, we introduce the scattering amplitude T a p t ys(ki, k 2 \ ks, k^) of triplet 

bosons where ki — (ki,k®). The ladder approximation satisfies a Bethe-Salpeter equation which is shown in Fig. [T] 
and written in Eq. (j22"|) . The scattering amplitude or self-energy for the two particle Green's function depends on 

the total energy and momentum of the incoming particles K = pi + p~2. The non- retarded and local character of U 
leads to T a p^s = ^ba-i&ps- The basic approximation made in the derivation of T(K) is that we neglect all anomalous 
scattering vertices, which are present in the theory due to existence of anomalous Green's functions. For the scattering 
amplitude shown in Fig. [TJ according to the Feynman rules in momentum space we can write (note p = (p , ~p) ) 

T a p, a p(piP2;P3Pi) = U(pi -p 3 ) 

+i(2tt)- 4 J [d 4 QU(Q~p 2 )G° aa (Q)G%( Pl +p 2 - Q)r a/S|a/S (p x +p 2 - Q,Q\pm)\ (22) 

In the above equation U is independent of momentum and energy. Consequently the T-function depends only on the 
sum of the incoming momentum and energy, and does not depend separately on the momentum and energy of the 
incoming particles. Therefore, p\ + p 2 = ps + p^ = K = (K , oj) which simplifies Eq. (|22p to 

T af3taP (K, uj) = U + i(2tt)- 4 J d 4 QUG° aa (Q)G° p p(K - Q)T aP . aP [[K, to). (23) 

However, the key observation is that all anomalous contributions are suppressed by an additional small parameter 
present in the theory - the density of the triplet excitation ni = ^ a (£„,-ia») — N~ l J2~tf a v q,a ~ 0T. Indeed, both 
terms of the anomalous scattering matrix are proportional to a which are neglected. By replacing the noninteracting 
normal Green's function (Eq. (fH)0 in the Bethe-Salpeter equation fEq. ([2"3")0 and taking the limit U — > oo we obtain 
the scattering matrix in the form (see Appendix A) 

iwj?, ")=-(4e <aU *- q * <y T p ) _1 - ( 24 ) 

q 

Now, we can calculate the single particle self-energy (Fig [5]) of bosons by utilizing the two particle self-energy (r) 
shown in Fig. [T] and obtained in Eq. ([24]^ . Because of the strong interaction between the triplet bosons we should 
carry out the expansion in Dyson's equation to infinite order. Therefore, S^ Q (fc) is written by 

/OO />OG 
d 4 P r af3nS (p,k;k,p)G° lf3 (p) + Y, d 4 P r a0n s(p,k;p,k)G° sp (p). (25) 

7 p -°° 05 J -°° 
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(Note that Y a p_~fS = Yap.apSa^Sps and U is frequency and momentum independent). For example, the x-component 
of the self-energy is written by 

E^(A) = 2(-) 4 / d±pT„(p + k)GUp) + (-r / d^pT^ip + k^ip) 

™ J-oo 27r J-oo 

/oo 
J 4 pT xy , xv (p + k)G yy (p). (26) 

In the above equation the first and third terms are similar. We integrate over the internal energy in complex plane 
on a contour in the upper half plane since G°(p) fEq. (|26p ) is anti-time ordered. Consequently, "E^ x (k, fco = u>) will be 
written by (N is the number of cells in the lattice) 

i . _ f°° -v 2 



/oo 
dp T XXtXX {p + k,p + k ] 

/oo 
dp T XZ}XZ {p + k,p + k 0/ 



,2 



V V 



In the dilute gas approximation there are other diagrams which are formally at most linear in n t (density of triplet 
bosons) but still numerically give contributions much smaller than Eq. I|27p . We should also consider the anomalous 
self-energy related to Hjj which is obtained from the vertex function in Eq. (|24|) . The anomalous self-energy (non- 
diagonal elements in the Eq. (fT7]) ) is written as (Fig [2]) 

S ^^E w ^ r (0,0). (28) 



The self-energy is proportional to ^2 q u q v q since the vertex function is independent of k. It is then proportional to 
the anomalous Green's function. It is negligible as compared with the normal self-energy (diagonal parts). 

We now consider the H3 contribution in the normal and anomalous parts of the self-energy. The normal part 
should be added to the self-energy due to Hjj. Since H 3 is much weaker than H\j it is sufficient to obtain the second 
order perturbation result of Dyson's series for each component of the normal Green's function. The formula for the 
self-energy contribution (either anomalous or normal) is quite lengthy and has been presented in Appendix B and the 
corresponding Feymann diagrams are shown in Fig. [3] The normal self-energy contribution of if 3 is proportional to 
u A which therefore dominates the anomalous one. 



V. EFFECT OF H 4 ON THE RENORMALIZATION OF THE SPECTRUM 



Because H4 is composed of quartic terms in the triplet operators its effect should be very small. It is therefore 
treated in mean field approximation by contracting the quartic operator products into all possible pairs. This is 
equivalent to take only the one loop diagrams (first order in J) into account. On mean field level we have 0\Oi = 
(Oi) O2 + (O2) 0\ — (O2) (0\) where each 0\ and O2 is a pair of boson triplet operator. We can write for each pair 
of operators 

(*U>) = /~ ^E eifc - ( ^ iMw0+G r(^) = ~£e*<*-*>t£ iaJ 

(*U>) = ^E eife ' ( ^" ii,) ^^' C 29 ) 

k 
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Thus, the effect of H4 is to renormalize A and B coefficients defined in H2 in the form 

A ktZ -^A ktZ +2J^-^2{v 2 g>x )^ q , 

1 

Bk,z > Bk :Z — 2J^fe— ^^(tt g ,x^q,g)£gi 

q 

Bk,(x,y) > B k ^ x , y )-J^ k —' S ^{5u q ^ Xi y)V q ,( x ,y)+U q ^V q ^)^ q . (30) 

1 

The renormalized coefficients fEq. (|3TJ)) ) will be considered to calculate the normal and anomalous self-energy 
which are independent of energy (nonretarded in time representation). The self-consistent solution of 
Eqs. (|24l27IC52IC53IC54IC55l30l2ip describes the quantum critical behaviour of this model which will be discussed in 
the following sections. 



VI. THE QUANTUM CRITICAL POINT AND THE GAP EXPONENT 

Close to the critical point quantum fluctuations exist over all length scales which define a scaling behavior for the 
physical quantities. The correlation length scales like £ ~ |Jj_ — J± c \ u where v is a critical exponent. This is related 
to the scaling behavior of the excitation gap in the Kondo singlet phase which vanishes like 

E g ~ |Jj_ - J^f, (31) 

as J± approaches its critical value J± c . Here cf> is called the gap exponent which is connected to the universality class 
of the quantum critical point. From general scaling arguments one expects (f> = vz where z is the dynamical critical 
exponent that determines the effective dimension -D e // of the model at T=0 according to -D e // = D + z. The spin 
excitation gap is defined by the energy of triplet excitations with x-polarization close to the antiferromagnetic wave 
vector qg. In its vicinity (\k — qo\ <C 1) the triplet dispersion can be approximated by 

LJ k<x = ^E? J+ 4(k-q ) 2 , (32) 

where c x is the spin- wave velocit y 13 i 19 ' 20 and qo = {tt,tt). There is no analytical expression for the spectrum of 
excitations, therefore we use numerical results to get the spin wave velocity. The slope of dispersion of the x- 
component excitations close to qo is calculated numerically which is the spin wave velocity c x . To find the energy gap 
we should consider the excitation energy at the wave vector qo 

Eg ~ Z q0!X (A q0tX - B 2 q X ), (33) 

where the renormalized constants have been obtained in the previous section 

Ak,x = ^(1 + A) + + S^(fc) + x {k) + ^LJ2( Z ^<* + Z «X,*)^ 

q 

B h . x = + E^(fc) + E^fc) - ^ ^2(Su q , x v g , x + u g , z v q>z )Z q . (34) 

q 

For the values of A q0iX and B q0:X at the critical point, A qoX = —B qoX holds. 

The energy gap in Eq. (|33[) vanishes at the quantum critical point ( Jj_ c ) and its behavior close to this point defines 
the scaling in Eq. (|3"T|) . We now look for the variation of the energy gap as J± deviates from J± c which is given by 
the variation of A q0jX and B qo ^ x with respect to J± deviation. The deviation of Jj_ from the critical point is defined 
by S S J± = J± — J± c . Therefore, close to critical point A qo x and B qox can be written as 

A qo ,x = A qo , x + l A^5 s J 1 _ + SE^tt, tt) + flE^Tr, tt) - ^ 

q 

Bq ,x = B^ x + dY? a>x {iz, tt) - ^ Yj /,.,<)„,,.,,;,., C r (35) 

q 
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where dX means the variation of X with respect to 8 S J±. If we substitute Eq. ([35|) into Eq. (|33|) and neglect terms 
quadratic in E g the variation of A qo and B qo must vanish 

dA qa = 1 A^ jL 5sJ± tt) + dV 3 n<x {iT,ir) - ^ Y.^- dv l^ = °. 

9S go = 9S^(7r, tt) - ^ ^;^,,^/,,,r, ( .,:^ = 0. (36) 

We now have to obtain the variation of each of the terms present in Eq. ([35| . In the first step, we calculate the 
variation of the self-energy related to Hjj which is written by 

j^^dv^ x T xx ^ x (q + q ,-iv q ^) + 3 J j—^v\ x dT XXiXX {q + q ,-u^ x ) 

;dvl z r xz , x 4q + q ,-uj qtZ ) + -^vl z dT xz , xz (q + q 0} -u) q , z ), (37) 



{2ir) 2 q ' z q J (2ir) 2 g 

where dv 2 z — 0. Indeed for < S < 1 the z-componcnt of the spectrum has a finite gap when the x-component 
becomes gapless at the quantum critical point. The main contribution to the first integral in Eq. (|37|) comes from the 
small momenta (q ~ E g <C 1) since 



and according to Eq. (|36|) the variation dA qx in this formula vanishes. We define the value of quantity X at the 
critical point by X c , Taking into account the first correction to the triplet density rib, the x-component of the vertex 
function can be written for small q (see Ref . ) 

pc 2 j^c 

F X x,xx(q, -u q , x ) » T c xx xx + XX '** c 2 °' X lug, (39) 

where T xx , xx (0, 0) = r£ x xx . The substitution of Eg. (f3"9"]l in Eq. (|3"T|) and replacing q — > S 9 /J in the first integral of 
Eq. (j37|) leads to the following equation for the variation of Hu self-energy 

o/lc Tp pc 2 AC Tp -i 

dJ%(„,*) = —^{T xx . xx + XX ZcC X J ] + f'^ nbS ^ +K x , xx n b S s J±. (40) 

In the second and fourth terms of Eq. (|3T|) we have T' a p Q/3 = STa i } - a j 3 j q '~^<') an( j rifj j g t ne density of triplet excitations 
at the critical point. Now, we consider the effect of H 3 on the gap exponent. According to Appendix B the variation 
of self-energy contribution of H 3 is given by 



on 



J 2 ( \ 1 + 2i> 2 



where x = —(4(5 + 2 + 2£ 2 ), (/? = —{28 + 1 + <5 2 ). Since in the vicinity of critical point 0J x ,q *C w 2 , ?0 the variation of 
Ea-pl]) gives 

-i ni2 ylc 

0( ) = g qo ' x , (42) 

+ W g> , 4(jB 2 + C 2 ((? _ 90 )2)1/2 K2 + ^2+^(^^)2)2 

The integral of Eq. (pt2"|) multiplied by 1 + 2w 2 a is proportional to E 2 In £" 9 which can be neglected compared with the 
first term in Eq. (|l0")) . Therefore, we can restrict to the variation of 1 + 2v 2 x . Then, the dominant contribution of 
Eq. (|4Tj) is given by 



A' 



JlB J 2 ((«o,»«o,*)x + (« 2 ) 2 + («S,») a )v) 
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We then calculate the last variation of the first expression in Eq.(|35 

5Z JA 



4 £(«W*£.& = ~^ E r ( 44 ) 
Let us define the following expressions 



o Ac yc 

° qo,x xx,xx 



4tt C 2 

^„/(K^)x + (K,) 2 + («S,*) 2 )v) 



CT = ^ r xz : xz n b$J± + T xxxx n b 5J x . (45) 
From the substitution of Eqs. (|44l43l40p into the first expression of Eq. (f3T>)) we get the following equation 

( + 2 3 + <t)5 s J± ^q ,x^xx,xx , S S J± 



Eb= x + e + v V- ^cKx + e + ^ -T)- (46) 

To find the gap exponent <j> we should consider E g — (£ s ./]_)• Indeed <j> is the smallest exponent that can be considered 
for Eg. Finally we obtain 

\AC pc 

i -| "^-qg.x-L xx,xx (A7\ 

9 ~ 47T C 2(A + + Ai)' 1 ' 

The last equation is obtained using uj x ^ qo <C For the isotropic case, Lu x ,q — ^>z,g = w goi where the expression 

for £^(7r, 7r) should be changed. The final result for the gap exponent of the isotropic case (6 = A = 1) is given by 

\A C „ r c 

J. - 1 » (48) 



4ttc 2 (\ + n + 6)' 



In the above equation we have 



AC _ Ac — AC 
^-q ,x ^q ,z — ^<Jo' 

A C n r c 
A = 90 



/it = 2 



9 ' 



^ J 2 (i6K^) + 8(K) 2 + K) 2 ) 

47rc 2 u;g 



(49) 



We have summarized the numerical values of the gap exponent for different anisotropics (S, A) in Tables Q] and HU 



VII. THE QUANTUM CRITICAL PHASE DIAGRAMS, NUMERICAL RESULTS 

Our approach is based on the strong coupling limit, Jj_/ J — > oo. In this limit the ground state has singlet character 
and a finite energy gap exists to the lowest excited triplet state. The increase of inter-site exchange coupling (J) 
or decrease of on-site exchange (J_i_) lowers the excitation gap which eventually vanishes. The position where the 
bosonic excitation gap vanishes defines the quantum critical point. At this point the condensation of a triplet takes 
place which induces the antiferromagnetic order. The single particle excitation energy should be obtained in a self 
consistent solution of Eqs. (|24l27IC52IC53IC54IC55l30l2ip . We should first replace 

Uk,a ► \J Zk, a Uk,a, Vk,a ► \fZ~k~aVk,a 
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in the self-energies of HjjjH^ and the renormalization expressions in Eg. ([30)) . From an initial guess for 
Zk a, S„ Q (fc, 0), Y* a a (k, 0) by using Eq. ([2"Tj) we obtain corrected excitation energy and the renormalized Bogoliubov 
coefficients (u, v). We repeat the procedure until the difference between the excitation energies in two consecutive 
steps is smaller than an acceptable error. The exponent of the gap is given by Eq. ([4"T[) which is calculated from the 
vertex function and spin wave velocity at the critical point and at the antiferromagnetic wave vector Q = (tt,tt). 
We will discuss the numerical results of our calculations in the next subsections for XY-case, i.e., 8 = and various 
sizes of on-site exchange anisotropy < A < 1 and likewise for A = 1 with various values of the inter-site exchange 
anisotropy < 5 < 1. 

A. XY-case: S = 

In the XY-case, the z-component of single particle excitation has approximately a dispersionless value ui z (k) = loq. 
For the other components of excitations (u> x = lu v ) the dispersion relation show a minimum at the reciprocal vector 
at the corner of the BZ, e.g., Q = (it, it). In Fig. ft] we have plotted the energy gap (E g /J±) versus control parameter 
Jj_/J. For all values of A the gap vanishes at the critical point (^j-) c , where the transition from Kondo singlet to 
the antiferromagnetic phase occurs. We have presented the numerical values of the critical point in Table U for 5 = 
and different values of intersite anisotropy, A = 0.0,0.2,0.4,0.6,0.8,0.9,1.0. In this case the critical points (J±/J) c 
have also been plotted in Fig. [6] versus A (solid line). In Table HI we have also compared the critical point with the 
mean field results^. The gap exponent for different anisotropies have also been presented in Table U and Fig. [7] Our 
results in this case show that the anisotropy in the local interaction (A) does not change the qualitative behavior of 
the Kondo necklace model on the 2D lattice. This is concluded from numerical values for the gap exponent ip which 
are independent of A. 

B. General anisotropic case 5^0 

To study the effect of a nonzero anisotropy 5 ^ in the itinerant part we consider only the isotropic case of the local 
Kondo interactions, i.e., J± x = J± z = J±, (A = 1). All of the three excitations show fc-dependence in this case. The 
minimum excitation energy at Q wave vector defines the energy gap. For Q < S < 1, we have Eg < Eg therefore we 
have plotted Eg /J± versus J±/J in Fig. [5j In contrast to the case of local anisotropy the effect of inter-site anisotropy 
on the critical point is weak and (Jj_/J) c changes only slightly. In Table [TT1 and Fig. [5] we have shown the calculated 
values of the critical point for different S. We have also presented the gap exponent for each anisotropy in Table HT1 and 
as a plot in Fig. [71 We observe a rapid decrease of <j> on approaching the isotropic case (6 = 1). In this case another 
soft mode (uj z ) is added to the excitations at the critical point and one needs a larger hopping strength (J) to reach 
the quantum critical point. Moreover, the symmetry changes from U(l) for < 8 < 1 to SU(2) at 8 = 1. Despite the 
weak dependence of the critical point on the anisotropy (<5) the critical exponent of the gap changes with 8 which is 
more pronounced at 8 = 1. However, the whole region of < <5 < 1 can be considered in a single universality class 
where the gap exponent changes slightly while the change of exponent at 5 = 1 signifies a different universality class 
by restoring the full spin rotational symmetry. 



TABLE I: The critical point (l 1 -^)^ at which the singlet gap vanishes for different values of intersite anisotropy (A) and 5 = 0. 
The second row shows the result from Green's function approach and the third row gives the mean field values for the critical 
poinli^. The gap exponent in the fourth row is obtained from the numerical evaluation of Eq.((47j. The accuracy of data is 
±0.005. 



A 


0.0 


0.2 


0.4 


0.6 


0.8 


0.9 


1.0 


(-j-) c (Green's function) 


3.01 


2.55 


2.17 


1.90 


1.72 


1.64 


1.55 


(^) c (mean field) 


2.86 


2.38 


2.04 


1.78 


1.59 


1.51 


1.43 


<f> (gap exponent) 


0.83 


0.83 


0.83 


0.82 


0.82 


0.82 


0.82 
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TABLE II: The critical point ((-j L ) c ) at which the singlet gap vanishes for different values of intersite anisotropy (S) and 
A = 1. The second row shows the result from Green's function approach and the third row is the mean field values for the 
critical point. The gap exponent in the fourth row comes from the numerical evaluation of Eqs. (|47l48p . The accuracy of data 
is ±0.005. 



s 


0.0 


0.2 


0.4 


0.6 


0.8 


0.9 


1.0 


(-j')c (Green's function) 


1.55 


1.54 


1.54 


1.52 


1.49 


1.47 


1.41 


(^j-)c (mean field) 


1.43 


1.43 


1.41 


1.39 


1.32 


1.30 


1.16 


4> (gap exponent) 


0.82 


0.82 


0.81 


0.81 


0.80 


0.78 


0.73 



VIII. DISCUSSION AND CONCLUSION 

In this work we have carried the analysis of quantum critical behavior of the 2D anisotropic Kondo necklace model 
beyond the previous mean field treatment. The constraint on the bosonic excitations has been implemented with 
the help of a hard core boson term at every site instead of applying a global constraint by introducing a chemical 
potential in the mean field approac h 14 ' 15 . 

The comparison of quantum critical points (J_i_/ J) c (<5 = l)as function of A in Table I shows that deviations of the 
two methods are quite small, up to 10% at the most. They are somewhat larger for the complementary case (A = 1) 
as function of 8 (up to ~ 17 %), especially when the isotropic case S = 1 is approached. The real difference and the 
advantage of the Green's function method appears when considering the critical exponent 4> of the excitation gap 
Eg. In the mean field treatment the exponent is always <p = 1 independent of the anisotropics. (Figs. 2 and 6 in 
Ref. [3). On the other hand the present Green's function approach clearly leads to nontrivial exponents (f> < 1 as 
may already be seen by the direct comparison of Figs. I4|5I with those of Ref. [3 mentioned above. The calculated 
critical exponents 4> are listed in Tables I and II and shown in Fig. [71 Generally they lie around </> ~ 0.80 — 0.83. 
The most remarkable feature is the rapid reduction of cj) when 6 approaches the isotropic point 6=1 where the 
universality class of the model changes from U(l)-xy to SU(2)-Heisenberg type. At this isotropic point we have 
<p ~ 0.73. Techniqually this means that three instead of two soft modes at the AF wave vector q appear which 
causes a rapid change in the gap exponent as described by Eqs. (I47I48[) . 

The critical value (J±/J) c and exponent <\> for the special isotropic case (5, A) = (1,1) has already been given 
in Ref. [l3|. From a numerical fit to the gap E g (Eq. I3T|) (J±/J) c = 1-39 and </> = 0.71 was obtained which agrees 
reasonably well with our values {J±/J) c — 1-41 and (f> = 0.73. In Ref. [H the isotropic case was also investigated 
for the true bilayer Hamiltonian. If we denote the inter-site coupling in the layer of localised spins (S^) by As then 
Ag = corresponds to the present KNM and As = J to the bilayer Hamiltonian. It was shown that the critical 
exponent cf> for the isotropic case does not depend on As- Applied to our anisotropic case we may conjecture that the 
critical exponents (/) given in Tables I and II will also be valid for the anisotropic bilayer model As = J although we 
have not done this calculation. 

The quantum Monte-Carlo results^ give v = 0.71, where v is the critical exponent of the correlation length close to 
quantum critical point of the 2D incomplete bilayer Heisenberg model (which is exactly the isotropic Kondo-necklace 
model). The dynamical exponent {£) relates the correlation length exponent {y) to the gap exponent (0) by <j) — zv. 
According to our results <f> = 0.73 which is very close to the result for v obtained in Ref. [23J, we conclude that the 
dynamical exponent of the 2D isotropic Kondo-necklace model is z = 1. As a consequence the effective dimension of 
the model at the QCP is D e ff = D + z = 3. This means that at the QCP the isotropic 2D KN model corresponds to 
the universality class of the 3D classical Heisenberg model. This was also suggested in Ref. [23|. 

It is also worthwile to discuss the different treatment of the bosonic constraint s]^ + ^ Q t\ a ti a = 1 (at every site 
i) in the two methods in more detail since it is the source of the difference in the position of the quantum critical 
point (J±/J) c mentioned above. Their essential difference has already been explained in Sect. [TTJ In the mean 
field approach a Lagrange term with a chemical potential /i incorporating the constraint is added to the quadratic 
Hamiltonian i?2- Then the singlet and triplet boson operators in the constraint are replaced by average amplitudes 
s and t a where the latter is identical to zero in the nonmagnetic phase. After diagonalization of H2 the singlet 
amplitude s and the chemical potential /1 are determined selfconsistently as function of J±/ J by minimizing the total 
energ y 14 ' 15 . This means that fi and s will be slowly varying functions of the control parameter J±/J. On the other 
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hand the present Green's function approach corresponds to fixing these two parameters to (Jj_/J) - independent 
constants given by s — 1 and fi = (Jj_/4)(2 + A) which correspond to the mean field values for J±/J = 0. Therefore 
on the level of H2 there is no constraint implemented in the present approach. It rather appears through adding 
the -ffy-term and computing its effect on the bosonic Greens functions in the limit U — > 00. One might speculate 
whether it would be an improvement to start from the self-consistent mean field solution for s and only then impose 
the further constraint with the Hy-term . However this cannot be justified easily. Once the constraints has been 
used in mean field level the fluctuations of the mean field singlet and triplet amplitudes should rather be unconstrained. 

However, the hard core repulsion on bosonic excitations which is imposed by Hjj, ensures that only one triplet can 
be excited on each bond. It justifies the dominant contribution to the excitation spectrum of the model. While, in the 
mean field approach there is no restriction to have more than one triplet excitation on each bond and the constraints 
s l s i + J2ocHa^ia = 1 should be statisfied only when averaged over all bonds. In other words, the Green's function 
approach may be expected to give an improved excitation spectrum. This is the reason why the critical exponent of 
the excitation gap obtained by Green's function approach is more accurate. 

The calculation of critical gap exponents presented here presents a contribution to understanding the quantum 
critical behavior of the anisotropic Kondo necklace model. An extension to the antiferromagnetic side of the quantum 
critical point would be desirable to allow for a full comparison with the mean field results in Ref. [3 However this 
will demand an even larger technical effort than in the present work. 



APPENDIX A 



The triplet scattering vertex in Eq.(|2"5|) is determined by the following equation 

T a ^ af5 (K,io) = - - . (B50) 

1 - iU(2n)-* J dZQdQ G° a (Q)G° pf3 (K - Q)T cM (K,cj) 

By substituting the normal noninteracting Green function Eq. (|14p in the integral in Eq. (|B50|) we have 

dQ G° aa {Q)G%{K - Q) = 

2 2 

I K -Q.fi K-QJ} s ,g 51 s 

Kq-Q - ujt^ p + id K a ~Q + ^_-^ <f3 - i5 

By using Cauchy's formula and extending the above integration over the complex plane Qo we compute the residues 
of the simple poles of the integrand. We found that the terms proportional to u?-> -» v\ and v\ -> u\ give zero 

K—Q,0 Q ,a K — Qfi Q ,a 

contribution because one pole is above the upper half plane and the other pole is below the lower half plane. After 
taking the limit U — ► oo then we obtain Eq. ([24|) for the scattering matrix. 
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APPENDIX B 



We have the following equations for x,z component of the normal and anomalous self-energies 
J 2 / 1 

*£,x(M) = ^T7 ( T [«9,*«,,x(Ufc +g ,, +t&H,, - 25^fc+, - ^ + 25^fc+ 9 ) 

g KUJ g-x + ^k+q.z 

+ (ul x U 2 k + qz + V 2 qx Vl + qz ){-5 2 il + q + SZgtk + g) 

+("fc+,,»«a,x + «fc+ g ,*«L)(-^fc + _ O 

+2Uk+ q , z Vk+q,zU q , x V q>x {^1 - 5£,k£,k+q - S,k£,q + £>£g£fc+g)] 

+ —~ L , Kv*(ut +gtX + vl +q . x mkt q 2&&+, - ^ + 2 5^*+«) 

+ K, 2 Wfe +9l2 + viz V k+ q ,z)(-^k+q + SZqtk+q) 
+Uk+ q ,xVk+q,x( u2 q,z + Vq,z)(£ktk+q ~ fk+q) 
+2Uk+q, x Vk+q, x U qiZ V qiZ (t,l - S£k£,k+q - £fc£g + £>£g£fc+g) 

+("fc+,,x«L + «fc+,,x«L)(-^ + - Ol) • (C52) 

For the anomalous x component self-energy we obtain 

J 2 ( 1 

S a,x(M) = ^T7^ (- — [Uq,xVq,x(Uk+q,z+Vk+q,z)(-^k£k+q+S 2 £k+q) 

+Mfe+g, z Wfc+g,«(M, ia; + ^ x )(25^fc^fe + g + ^ ~ £fc£g - 2S£ q £ k+q ) 

+ («fc+g,z«g,x + "fc+9,z«g,x)(-^fc + - <^g£fc+g + £fc£g) 

+ T" [ U q,zVq,z( U l+q,x + V l+q,x)(-^k+q + Cfe+g) 

+u fe+ g, a; w fe+9 , 2: (u^ + «g,J(266+g + £ 2 £g - <5£fc£g - 2^g6+g) 

+ 2Mfe +9>a: Wfe +g>a: Ug iZ Wg i2 (^ +9 - ft^fc+g + £fc + (5 2 ^ - 2<5£g£fc) 

+ ("fc+,,x"fl,z + W fc+g,x W g, z)(-^fc + &&+9 - ^gCfe+g + ^fe^g)]) • (C53) 

Similarly we can write the following expression for z component of the normal self-energy 

K,z(k,0) = — ^ ,, — [Uq,xVq,x(Uk+q,x + 4+q,x)(^q ~ 25&&+, ~ £g + 2£g&+g) 

g 1- x ' w fe+9,x 

+ (^g,x U fc+g,x + ^g,x^+g,x)(-^fe+g + 6?6+g) 

+ («fe+g,x«g,x + U 2 k+q , x vl x )(-S 2 Z 2 k + - Cg) 

+U k+q , x V k+q . x {u 2 q x + ^J^fc^+g - <5 2 Cfe+g) 

+ 2u k+q ^ x V k+qiX U q . x V q/x (S 2 ( 2 k - S^k^k+q - S£ k £,q + £g£fc+g)]- (C54) 

For the corresponding anomalous self-energy we will have 

^L(fc,0) = "Tr^" [Uq.xV q , x {ut + q , x + Vl +qix ){-5^k+q +il+q) 

q 1- x UJ k+ q ,x 
+U k +q,xVk+q.x{ u %x + v lx)( 2 ^k+q + £ q ~ 5£ k £q ~ 2£g£fc+g) 
+ ( V k+q,x u2 q,x + U l+q,x V q,x)(-<i 2 t.k + ^k£,k+q ~ £,qtk+q + 5£ k £ q ) 

+2u k+q!X V k + q ,xU q ,xV q ,x(tl +q ~ iqik+q + S 2 £,l + (% - 26£ q £ k )]- (C55) 
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FIG. 1: Ladder diagram for the triplet boson scattering amplitude T a p^s (k\, k^\ kz, ki). 




FIG. 2: Diagrams for the normal and anomalous single particle self-energy due to hard core term Hu and using the two-particle 
scattering amplitude in Fig. [1] 

E 3 - 

n 

^ a 

FIG. 3: One loop diagrams for normal and anomalous self-energy arising from the three point interaction H3. 
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0.2 0.4 0.6 0.8 
A(solid) or 5(dash) 

FIG. 6: Dependence of the critical point (J±/J) c on the local (A, solid line) and itinerant (8, dashed line) anisotropies. For 
the solid line 5 = and for the dashed line A = 1. The numerical values are given in Tables I and II. 
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FIG. 7: Dependence of the critical exponent (<f>) on the local (A, solid line) and itinerant (S, dashed line) anisotropies. For the 
solid line 8 = and for the dashed line A = 1. The numerical values are given in Tables I and II. 



